Dominance

This plot shows the dominance value of each genotype (given by within individual allele frequency) as a function of the dominance scalar. The purpose was to write a function that would take a single dominance value that could be specified for any ploidy and return a comparable set of dominance coefficients for heterozygous genotypes.

getDomCoefs = function(dominance, ploidy = 100){
  if (dominance == 0){
    coeffs = rep(0, ploidy - 1)
  }
  else if (dominance == 1){
    coeffs = rep(1, ploidy - 1)
  }
  else {
    dom = (1 - dominance) - 0.5
    dom = abs((dom * 10 + sign(dom)) ^ sign(dom))
    coeffs = (seq(1, ploidy - 1) / ploidy) ^ dom
  }
  return(coeffs)
}

doms = c(seq(0,1, 0.1))
Doms = as.data.frame(sapply(doms, getDomCoefs))
Doms = rbind(c(0,0,0,0,0,0,0), Doms, c(1,1,1,1,1,1,1))
colnames(Doms) = doms
Doms = cbind(Doms, data.frame("i/k"=seq(0,100) / 100))
mDoms = melt(Doms, id.var = "i.k")
ggplot(mDoms, aes(x = i.k, y = value, color = variable)) + geom_line() + scale_color_discrete(name = "Dominance\nScalar") + theme_bw() + xlab("i / k") + ylab("Dominance Coefficient")

# Example

# dominance=c(seq(0.1,0.9,0.1))
# > dominance
# [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

# Centers around 0.5 (additivity)
# dom = (1 - dominance) - 0.5
# > dom
# [1]  0.4  0.3  0.2  0.1  0.0 -0.1 -0.2 -0.3 -0.4

# dom1 = abs((dom * 10 + sign(dom)) ^ sign(dom))
# > dom1
# [1] 5.0000000 4.0000000 3.0000000 2.0000000 1.0000000 0.5000000 0.3333333 0.2500000 0.2000000

# > getDomCoefs(dominance = 0.2, ploidy = 4)
# [1] 0.00390625 0.06250000 0.31640625
# > getDomCoefs(dominance = 0.5, ploidy = 4)
# [1] 0.25 0.50 0.75

Sweep trajectories

These plots show a primary motivating result for this study, in that higher ploidies show increased fixation times of sweeping beneficial alleles regardless of the dominance scenario being evaluated.

#Function creates trajectories
getTraj = function(s, dom, ploidy, start_freq, N = -9, end_freq = 0.99, maxGens = 9999, maxTries = 10){
  dom_coeffs = getDomCoefs(dom, ploidy)
  attempts = 0
  p = start_freq
  k = ploidy
  gen = 0
  while (attempts <= maxTries & p < end_freq & gen <= maxGens){
    df = data.frame("s" = s, "dom" = dom, "gen" = 0, "freq" = start_freq, "dp1" = 0, "w.bar" = 0, "ploidy" = ploidy)
    p = start_freq
    dp1 = 0
    gen = 0
    het_fits = 1 + dom_coeffs * s
    fits = c(1, het_fits, 1 + s) / (1 + s)
    while (p < end_freq & p > 0 & gen <= maxGens){
      q = 1 - p
      i = seq(0, k)
      Num = sum(choose(k, i) * ((k - i) / k) * (p ^ (k - i)) * (q ^ i) * rev(fits))
      Den = sum(choose(k, i) * (p ^ (k - i)) * (q ^ i) * rev(fits))
      p_prime = Num / Den
      if (N != -9){
        p_prime = sum(rbinom(N, k, p_prime)) / (k * N) 
      }
      df = rbind(df, c(s, dom, gen, p, dp1, Den, k))
      dp1 = p_prime - p
      p = p_prime
      gen = gen + 1
    }
    df = rbind(df, c(s, dom, gen, p, dp1, Den, k))
    attempts = attempts + 1
  }
  if (attempts > maxTries & nrow(df) <= maxGens){
    print(paste("Beneficial allele was lost due to drift for", maxTries, "consecutive attempts"))
  }
  if (nrow(df) > maxGens){
    print(paste("Beneficial allele did not fix before the maximum number of generations (", maxGens, ")."))
  }
  return(df[-1,])
}

# Param set to create trajectories for
params = data.frame("s" = rep(0.1,9), "dom" = c(rep(0.5,3), rep(0,3), rep(1,3)), "ploidy" = rep(c(2,4,8),3), "start" = rep(0.05, 9))

# Create trajectories
dat1 = foreach(i = 1:nrow(params), .combine=rbind) %do% {
  getTraj(params[i,]$s, params[i,]$dom, params[i,]$ploidy, params[i,]$start)
}

## [1] "Beneficial allele did not fix before the maximum number of generations ( 9999 )."
## [1] "Beneficial allele did not fix before the maximum number of generations ( 9999 )."
## [1] "Beneficial allele did not fix before the maximum number of generations ( 9999 )."
## [1] "Beneficial allele did not fix before the maximum number of generations ( 9999 )."

#plot
ggplot(dat1, aes(x=gen, y=freq, linetype=as.factor(dom), color=as.factor(ploidy))) + geom_line() + scale_color_discrete(name="Ploidy") + scale_linetype_discrete(name="Dominance Scalar", labels=c("0.1", "0.5", "0.9")) + xlim(0,300) + theme_bw() + xlab("Generation") + ylab("Allele Frequency")

## Warning: Removed 40024 rows containing missing values (geom_path).

Fitness

First plot shows reason why allele frequency change is faster in lower ploidies; increased variance in fitness at any given frequency. The second plot just shows mean fitness as a function of allele frequency and is not terribly interesting.

getFits = function(s, dom, ploidy){
  dom_coeffs = getDomCoefs(dom, ploidy)
  k = ploidy
  dat = foreach(i = 1:51, .combine=rbind) %do% {
    p = seq(0,1,1/50)[i]
    q = 1 - p
    het_fits = 1 + dom_coeffs * s
    fits = c(1, het_fits, 1 + s) / (1 + s)
    i = seq(0, k)
    freqs = choose(k, i) * (p ^ (k - i)) * (q ^ i)
    w.bar = sum(freqs * rev(fits))
    var.w = sum((freqs * (rev(fits) - w.bar) ^ 2) / length(freqs))
    data.frame('s' = s, 'dom' = dom, 'ploidy' = ploidy, 'freq' = p, 'w.bar' = w.bar, 'var.w' = var.w)
  }
  return(dat)
}

dat2 = foreach(i = 1:nrow(params), .combine=rbind) %do% {
  getFits(params[i,]$s, params[i,]$dom, params[i,]$ploidy)
}

ggplot(dat2, aes(x=freq, y=var.w, color = as.factor(ploidy), linetype=as.factor(dom))) + geom_line() + scale_color_discrete(name="Ploidy") + scale_linetype_discrete(name="Dominance", labels=c("0.1","0.5","0.9")) + xlab("Allele Frequency") + ylab("Variance in Fitness") + theme_bw()

ggplot(dat2, aes(x=freq, y=w.bar, color = as.factor(ploidy), linetype=as.factor(dom))) + geom_line() + scale_color_discrete(name="Ploidy") + scale_linetype_discrete(name="Dominance", labels=c("0.1","0.5","0.9")) + xlab("Allele Frequency") + ylab("Mean Fitness") + theme_bw()

Trajectory effect

Here, we are keeping fixed the ploidy specified for the coalescent simulations (fixed at 8 in this example), and we are varying the ploidy under which the sweep trajectories are simulated. For some parameter sets, higher ploidies and thus slower sweeps results in a visible difference in terms of the reduction of diversity at and around the beneficial allele. The effect is more pronounced for weaker selection and for partial dominance/recessivity. The breadth of the dips is also notably lesser for higher ploidies. Thus, if chromosome number in a population is held constant, higher ploidies exhibit reduced signals associated with selective sweeps.

ms2 %>% filter(dom %in% c("0.1","0.9","0.5") & N == 10000 & fuseGen == 1 & sampGen ==1 & mu==1e-08 & sim_ploidy == 8) %>% ggplot(., aes(y=Pi.1 * 100, color = as.factor(traj_ploidy), x = mid/100000)) + geom_smooth(method="loess", size =0.5, span=0.1, se=F) + scale_x_continuous(breaks=c(1,5,9)) + xlab("Position (100 kb)") + ylab("Diversity (%)") + scale_color_discrete(name="Trajectory\nPloidy") + facet_grid(vars(s, recomb), vars(dom),labeller = label_both, scales = "free_y") + theme_bw()

Factor-of-2 effect

This is the complement to the previous figure. We use the same sweep trajectory, but we vary the ploidy specified for the coalescent simulation.

ms2 %>% filter(dom %in% c("0.1","0.9","0.5") & fuseGen == 1 & sampGen ==1 & mu==1e-08 & N == 10000 & traj_ploidy == 4) %>% ggplot(., aes(y=Pi.1 * 100, color = as.factor(sim_ploidy), x = mid/100000)) + geom_smooth(method="loess", size =0.5, span=0.1, se=F) + scale_x_continuous(breaks=c(1,5,9)) + scale_y_continuous(breaks=c(0.05,0.15)) + xlab("Position (100 kb)") + ylab("Diversity (%)") + scale_color_discrete(name="Simulation\nPloidy") + facet_grid(vars(recomb,s), vars(dom),labeller = label_both) + theme_bw()

Combined effects

Higher ploidies result in deeper and broader dips in diversity relative to diploids. The effect of recombination rate is not striking.

ms2 %>% filter(dom %in% c("0.1","0.9","0.5") & fuseGen == 1 & sampGen ==1 & mu==1e-08 & N == 10000 & sim_ploidy == traj_ploidy) %>% ggplot(., aes(y=Pi.1 * 100, color = as.factor(traj_ploidy), x = mid/100000)) + geom_smooth(method="loess", size =0.5, span=0.1, se=F) + scale_x_continuous(breaks=c(1,5,9)) + scale_y_continuous(breaks=c(0.05,0.15)) + xlab("Position (100 kb)") + ylab("Diversity (%)") + scale_color_discrete(name="Ploidy") + facet_grid(vars(recomb,s), vars(dom),labeller = label_both) + theme_bw()

Tajima’s D

Polyploids recover new mutations more quickly, but they take longer to drift to higher frequencies. ThetaW goes up but ThetaPi lags…

# ms2 %>% filter(dom %in% c("0.1","0.9","0.5") & s==0.01 & fuseGen == 1 & mu==1e-08 & N == 10000 & sim_ploidy == traj_ploidy & recomb == 1e-08) %>% ggplot(., aes(y=Tajima.D_1, color = as.factor(sim_ploidy), x = mid/100000)) + geom_smooth(method="loess", size =0.5, span=0.1, se=F) + scale_x_continuous(breaks=c(1,5,9)) + xlab("Position (100 kb)") + ylab("Tajima's D") + scale_color_discrete(name="Ploidy") + facet_grid(vars(sampGen), vars(dom),labeller = label_both) + theme_bw()

ms2 %>% filter(dom == "0.5" & s==0.01 & fuseGen == 1 & mu==1e-08 & N == 10000 & recomb == 1e-08) %>% ggplot(., aes(y=Tajima.D_1, color = as.factor(sim_ploidy), x = mid/100000)) + geom_smooth(method="loess", size =0.5, span=0.1, se=F) + scale_x_continuous(breaks=c(1,5,9)) + xlab("Position (100 kb)") + ylab("Tajima's D") + scale_color_discrete(name="Ploidy") + facet_grid(vars(sampGen), vars(traj_ploidy),labeller = label_both) + theme_bw()

## Warning: Removed 204 rows containing non-finite values (stat_smooth).

Fst

Neutral divergence in allele frequencies increases Fst more quickly in lower ploidies, thus erasing ‘outlier’ signal due to sweep

ms2 %>% filter(dom %in% c("0.1","0.9","0.5") & s==0.01 & fuseGen == 1 & mu==1e-08 & N == 10000 & sim_ploidy == traj_ploidy & recomb == 1e-08) %>% ggplot(., aes(y=fst, color = as.factor(sim_ploidy), x = mid/100000)) + geom_smooth(method="loess", size =0.5, span=0.1, se=F) + scale_x_continuous(breaks=c(1,5,9)) + xlab("Position (100 kb)") + ylab("Fst") + scale_color_discrete(name="Ploidy") + facet_grid(vars(sampGen), vars(dom),labeller = label_both) + theme_bw()